Draft version August 24, 201 1 

Preprint typeset using ET^X style emulateapj v. 1 1/10/09 



A NEW STELLAR MIXING PROCESS OPERATING BELOW SHELL CONVECTION ZONES FOLLOWING 

OFF-CENTER IGNITION 

M. Mocak 1 , Casey A. Meakin 2 ' 3 , E. Muller 4 & L. Siess 1 
Draft version August 24, 2011 

ABSTRACT 

During most stages of stellar evolution the nuclear burning of lighter to heavier elements results in a ra- 
dial composition profile which is stabilizing against buoyant acceleration, with light material residing above 
heavier material. However, under some circumstances, such as off-center ignition, the composition profile 
resulting from nuclear burning can be destabilizing, and characterized by an outwardly increasing mean molec- 
ular weight. The potential for instabilities under these circumstances, and the consequences that they may have 
on stellar structural evolution, remain largely unexplored. In this paper we study the development and evolution 
of instabilities associated with unstable composition gradients in regions which are initially stable according 
to linear Schwarzschild and Ledoux criteria. In particular, we explore the mixing taking place under various 
conditions with multi-dimensional hydrodynamic convection models based on stellar evolutionary calculations 
of the core helium flash in a 1.25 M star, the core carbon flash in a 9.3 M star, and of oxygen shell burning 
in a star with a mass of 23 M© . The results of our simulations reveal a mixing process associated with regions 
having outwardly increasing mean molecular weight that reside below convection zones. The mixing is not due 
to overshooting from the convection zone, nor is it due directly to thermohaline mixing which operates on a 
timescale several orders of magnitude larger than the simulated flows. Instead, the mixing appears to be due to 
the presence of a wave field induced in the stable layers residing beneath the convection zone which enhances 
the mixing rate by many orders of magnitude and allows a thermohaline type mixing process to operate on a 
dynamical, rather than thermal, timescale. The mixing manifests itself in the form of overdense and cold blob- 
like structures originating from density fluctuations at the lower boundary of convective shell and "shooting" 
down into the core. They are enriched with nuclearly processed material, hence leaving behind traces of higher 
mean molecular weight. In these regions we find that initially smooth composition gradients steepen into stair- 
step like profiles in which homogeneous, mixed regions are separated by composition jumps. These step like 
profiles are then seen to evolve by a process of interface migration driven by turbulent entrainment. We discuss 
our results in terms of related laboratory phenomena and associated theoretical developments. We also discuss 
the degree to which the simulated mixing rates depend on the numerical resolution, and what future steps can 
be taken to capture the mixing rates accurately. 

Subject headings: Stars: evolution - hydrodynamics - convection - mixing 



1. INTRODUCTION 

The hydrodynamic approach to model certain phases of 
stellar evolution, like e.g., a nuclear core flash, by numeri- 
cally solving the Euler or Navier-Stokes equations is essen- 
tially built upon first principles in physics, and thus (almost) 
parameter-free. This approach is advantageous compared to 
(ID) stellar evolutionary calculations when modeling phe- 
nomena related to buoyancy or shear driven flow instabilities 
which are inherently multidimensional in nature. A shortcom- 
ing of this approach is that it remains impossible to simulate a 
star over evolutionary timescales. Instead, multi-dimensional 
hydrodynamic calculations provide insight into the nature of 
mixing processes operating over fluid dynamical timescales 
which can be used to develop more comprehensive mixing 
models to be used in evolutionary codes. The classic mixing 
length theory (MLT) for thermal convection, which remains 
the standard treatment of convection in stellar evolution cal- 
culations, is an example of one such mixing model based on, 
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in this case, a very crude picture of turbulent flow. A grow- 
ing body of work presenting simulations of stellar interiors 
(e7g. [D earborn et a l.|2006[|H erwig et al.|2006|[2QTT| [M eakin 
& Arnett 2007 j |Arnett et al.|2009[ |Mocak et al.|2U0gr f2009 ■ 
|2010| ) offers the promise of moving beyond crude models like 
MLT by providing theorists with detailed information about 
the full non-linear development of fluid instabilities under re- 
alistic astrophysical conditions. 

We build on this work in the present paper by providing a 
detailed account of a new mixing process that we have dis- 
covered in simulations of shell convection. In particular, we 
present calculations of both helium and carbon shell convec- 
tion following off-center ignition. In both of these evolu- 
tionary phases we identify a coupling between the convec- 
tive shell and the layers which reside below the shell. This 
coupling leads to mixing below the burning shell that has po- 
tentially significant consequences for the star's structure and 
evolution. 

The paper is organized as follows. In Sect.[2]we described 
our hydrodynamics code and input data and present some gen- 
eral features of the simulated flows. In Sect.[3]we provide a 
detailed account of the new mixing process that we have iden- 
tified and discuss its properties in the context of other pro- 
cesses commonly encountered in stellar interiors, including 
penetration and overshoot, semiconvection, and salt finger- 
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TABLE l 

Some characteristic properties of our initial models: total stellar mass 
m, stellar population, metal content z, mass of the mapped model m m , 
outer (inner) radius r m of the mapped model, and nuclear energy 
production rate l m of the mapped model, respectively. 



TABLE 2 

Some properties of the 2D and 3D simulations: number of grid points in r 

(N r ), 6 (Ng), AND (N(p) DIRECTION, RADIAL GRID RESOLUTION Ar, ANGULAR GRID 
RESOLUTION AO AND A0 IN thetd AND p/zj'-DIRECTION, AND MAXIMUM 
EVOLUTIONARY TIME t max OF THE SIMULATION, RESPECTIVELY. 
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ing. Contact is made with relevant laboratory, geophysical, 
and fluid mechanics literature. We summarize our findings in 
Sect.|4j and indicate future research directions. 

2. INITIAL DATA AND SIMULATION PROPERTIES 
2.1. Initial data 

We have simulated the reactive hydrodynamic evolution 
for three phases of stellar evolution. The initial data used 
for our calculations include: (1) the helium core of a 1.25 
M Q star during the peak (maximum core luminosity) of the 
core helium flash comp uted with the GARSTEC code |Weiss 
|& Schlatfl|2000| 12007] ), (2) the carbon-oxygen (C-O) core of 
a 9.3 M G star at peak o f the core ca rbon flash computed with 
the STAREVOL code ( [Siessl[20Q6] ), and (3) the core of a 23 
M© star during oxygen shell burning (Meakin & Arnett 2007 ) 
computed with the TYCHO code ( |Arnett et al.||2010| ). The 
thermodynamic structure of all three initial models is qualita- 
tively similar exhibiting an off-center temperature maximum 
due to nuclear burning in a partially electron degenerate stel- 
lar core (see Fig.[T]). Some additional properties of the stellar 
models are included in Table [T] 

The core helium flash occurs after central hydrogen exhaus- 
tion in the semi-degenerate helium core of low-mass stars 
(0.7 M© < M < 2.2 M© ) due to an off-center ignition of he- 
lium by the triple-a. The core carbon flash ensues after central 
helium exhaustion and is also characterized by off-center car- 
bon ignition in a semi-degenerate carbon-oxygen (CO) core 
of a rather massive star (7M < M < 11M ), the domi- 
nant nuclear reactions being 12 C ( 12 C, a) 20 Ne and 12 C( 12 C, 
p) 23 Na followed by 16 0(a, y) 20 Ne. Oxygen shell burning, 
which is typical for massive stars (M > 12 M© ) close to core 
collapse, is not ignited off-center, but follows an epoch of 
core oxygen burning which leaves behind a silicon- sulfur rich 
semi-degenerate core. Because of the steep temperature de- 
pendence of the nuclear energy production rates, convection 
develops in the burning shell. 

During the core helium flash the convective shell is en- 
riched mainly in 12 C which results in a negative mean molec- 
ular weight gradient below its base (Fig.[l]a) |^] The situa- 
tion is similar during the core carbon flash, where the nuclear 
ash from carbon burning ( 20 Ne and 16 O mainly) increases the 
mean molecular weight inside its convective shell relative to 
that of the unburned inner CO core (Fig.[T]b). Oxygen shell 
burning on the other hand, which follows oxygen core burn- 
ing, results in a positive composition gradient, i.e., the mean 
molecular weight increases in the direction of gravity at the 
bottom of the convective shell (Fig.[T]c). 

2.2. Properties of hydrodynamic simulations 



5 The gradient grows for roughly 10000 years since the onset of the core 
helium flash. 
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The numerical simulations were performed with a modi- 
fied v ersion of the hydrodynamic code Herakles ( Mocak et al. 
2008) . The code employs the PPM reconstruction scheme 
( |Cole lla & Woodward 1984) a Rie mann solver for real gases 
according to Colella & G laz (198 4)), and the consist ent multi- 
fluid advection scheme of [Plewa & Mtiller| ( |1999| ). Self- 
gravity, thermal transport and nuclear burning are included in 
the code. Nuclear reaction networks are generated using the 
REACLIB library developed by Thielemann (private commu- 
nication). In Table[2] we summarize some characteristic pa- 
rameters of our 2D and 3D hydrodynamic simulations. 

The core helium flash simulations were performed on an 
equidistant spherical polar grid in 2D (r, 0) and 3D (r, 0, 0). 
While the angular grid of the axisymmetric 2D simula- 
tions hefl.2d.2 and hefl.2d.3 covered the full angular range, 
i.e., 0° < < 180°, the 2D simulation hefl.2d.l was restricted 
to the angular region 30° < < 150°. The 3D simula- 
tion hefl.3d was performed within an angular wedge given 
by 30° < 6 < 150° and -60° < (p < +60°, respectively. 
This allowed us to simulate the 3D model with a reasonable 
amount of computational time using a radial and angular res- 
olution comparable to that of the 2D model hefl.2d.l for sev- 
eral convective turnover timescales. In all core helium flash 
simulations the radial grid ranged from r = 2 x 10 8 cm to 
r - 1.2 x 10 9 cm. Boundary conditions were reflective in all 
directions except for simulations hefl.2d.l and hefl.3d, where 
we utilized periodic boundary conditions in angular direc- 
tion^) to avoid a numerical bias in case of wide angular con- 
vective structures. Abundance changes due to nuclear burning 
during the He-flash are described by a reaction network con- 
sisting of the four a-nuclei 4 He, 12 C, 16 O, and 20 Ne coupled 
by seven reactions. 

The core carbon flash simulation cafl.2d was performed on 
a 2D equidistant spherical polar grid (r, 6) covering a 90° an- 
gular wedge centered at the equator (6 = 90°), and a radial 
grid ranging from r = 3 x 10 8 cm to r = 1 x 10 9 cm. Boundary 
conditions were reflective in radial direction and periodic in 
angular direction. Abundance changes due to nuclear burn- 
ing were described by a reaction network consisting of 1 H, 
4 He, 12 C, 14 N, 16 O, 20 Ne, 22 Ne, 23 Na, and 24 Mg coupled by 17 
reactions. 

The oxygen shell burning simulation oxfl.2d was performed 
on a 2D equidistant spherical polar grid covering a 90° wedge 
centered at the equator, and a radial grid ranging from r = 
2 x 10 8 cm to r = 1 x 10 9 cm. Boundary conditions were 
reflective in radial direction and periodic in angular direction. 
Abundance changes due to nuclear burning were described by 
a reaction network consisting of neutrons, 1 H, 4 He, 12 C, 16 O, 
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Fig. 1. — Temperature T (solid) and mean molecular weight // (dotted) as a function of radius for the initial core helium flash model (a), the core carbon flash 
model (b), and the oxygen shell burning model (c), respectively . In each of these panels the region of interest below the base of the convection zone is highlighted 
by the shaded vertical strip. 
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Fig. 2. — Snapshots of the relative angular fluctuations of the mass fraction of 4 He during the core helium flash for model hefl.2d.3 at t ~ 12000 s (a), of 12 C 
during the core carbon flash for model cafl.2d at t ~ 682 s (b), and of 16 during oxygen shell burning for model oxfl.2d at t ~ 940 s (c), respectively. The 
fluctuations are defined by AX( A N) = 100 x [X( A N) - (X( A N)) e ] I (X( A N)) e , where A N e { 4 He, 12 C, 16 0}, and {) e denotes the angular average at a given 
radius. 
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Ne, 23 Na, 24 Mg, 28 Si, 31 P, 32 S, 34 S, 35 C1, 36 Ar, 38 Ar, 39 K, 



40 Ca, 42 Ca, 44 Ti, 46 Ti, 48 Cr, 50 Cr, 5Z Fe, M Fe, and 5b Ni coupled 
by 76 reactions. 

The global character of the flow is illustrated in Figure [2] 
which shows a snapshot of the composition fluctuations in the 
developed convective flows from the hefl.2d.3 (helium burn- 
ing), cafld.2d (carbon burning), and oxfl.2d (oxygen burning) 
calculations. The basic topology of the flow is that of a tur- 
bulent convective shell surrounded by stably stratified layers 
which host internal waves. In these 2D simulations the shell 
convection consists of interacting plumes and convective rolls 
which span the depth of the shell. In 3D the flow is better 
described by plumes than rolls, since vortices are no longer 
pinned to mer idional planes by the ge ometry of the calcula- 



tion (see e.g. |Meakin & Ar nett 2007). In both the 2D and 



the 3D simulations, the internal waves can be seen as nar- 
row (in radius) bands extending in angle above and below the 
turbulent convective shell. The lack of a strong signature of 



internal waves above the convective shell in the helium and 
oxygen burning models in Figure [2] is a result of the compo- 
sition profile lacking a gradient there. The signature of waves 
is present in other fields which have gradients, such as the 
pressure, density, and velocity fields. 

Helium Flash Convection. — The hydrodynamic properties 
of shell con vection during t he cor e helium flash are described 
in detail in |Mocak et ak] fl2008| |2009| ). Convection starts 
early (t < 1000 s) and quickly extends over the whole con- 
vectively unstable region. In axisymmetry (i.e., 2D), the 
convection is characterized by fast and large circular vor- 
tices, while in 3D the convective flow is less ordered showing 
slower and smaller turbulent features. The convective veloci- 
ties in our 3D model match those predicted by MLT quite well 
(|v| ~ v M lt ^ 1 X 10 6 cm s" 1 ), whereas the velocities in our 
2D models exceed those by up to a factor of four. Turbulent 
entrainment leads to a growth of the width of the convection 
zone on a dynamic timescale in both the 2D and 3D models 
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Fig. 3. — Maps of the relative angular fluctuation in carbon mass fraction AX( U C) = 100 x (X( U C) - (X( u C)) e ) / (X( u C)) e at the bottom of the He-flash 
convection zone, where ()g denotes the horizontal average at a given radius. From left to right, the following models are shown: (left) hefl.2d.3 at t = 63430 s, 
(middle) hefl.2d.l at t = 6000 s , and (right) a meridional plane of the 3D model hefl.3d at t = 6000 s. The vertical solid line marks the bottom boundary of the 
convection zone which is equal to the position of T max , and the dashed line gives the location from where our mixing begins. 



due to an exchange between the potential energy contained in 
the stratified layers at the boundaries of the convection zone 
and the kinetic energy of the turbulent flow inside the convec- 
tion zone. 

Carbon Flash Convection. — An analysis of the hydrody- 
namic properties of shell convection during the core carbon 
flash, based on our 2D model carl. 2d, shows that the con- 
vective flow is dominated by small circular vortices. The an- 
gular averaged amplitudes of velocities inside the convection 
zone |v| ~ 4 x 10 6 cm s" 1 exceed those predicted by MLT 
Vmlt ~ 1.5 x 10 5 cm s" 1 by about an order of magnitude. 
From authors experience, corresponding 3D convection (not 
calculated) would be characterized by smaller velocities by a 
factor from 2 to 5. The operation of turbulent entrainment at 
the convective boundaries has also been found but not inves- 
tigated in detail in this paper. 

Oxygen Shell Burning. — The hydrodynamic properties of 
shell conve ction during oxygen she ll burning are discussed 
in detail in |Meakin & Arnett| ( |2007] ). Our 2D model oxfl.2d 
reproduces approximatelly the results of these authors. The 
angular averaged amplitudes of the convective velocities in- 
ferred from model oxfl.2d are |v| ~ 1 x 10 7 cm s" 1 , which 
is roughly equivalent to the velocity v M lt predicted by MLT. 
Turbulent entrainment has been detected in model oxfld.2d, 
too, but we hav e not further analyzed th is phenomenon; for 
more details see Meakin ^& Arnett| ( |2007| ) . 

3. RESULTS 
3.1. General Features of the Mixing 

In this paper we focus our analysis on a mixing process 
which operates within the stably stratified layers residing be- 
neath the shell convection zones in both the carbon- and 
helium-flash simulations. This mixing process is not observed 
to operate below the shell convection zone in the oxygen burn- 
ing case. In this subsection we describe the salient features of 
the mixing for the two shell flash calculations, and contrast 
these with the oxygen shell burning case. 

Helium Flash Convection. — In the helium-flash convec- 
tion simulations we find a mixing process taking place be- 
low the temperature maximum immediately after convection 
starts. The mixing manifests as the formation of blob-like 
structures in the stable layer just below the temperature maxi- 
mum, which fall inward toward the stellar center. These struc- 
tures are present in both the 2D and 3D simulations and can 
be seen in the snapshots presented in Figure [3] which shows 



the fluctuations in the helium burning product C 12 for three of 
our He-flash simulations. 

The layers just beneath the base of the convection zone 
show the presence of overdense sinking blob-like structures 
which are enriched with higher mean molecular weight \i 
material than the ambient matter into which they penetrate, 
thereby creating a compositional step as they redistribute the 
yu|^] The beginning of new compositional step essentially 
coincides with the layer from which the mixing launches. 
Whereas it almost overlaps with the bottom of convection 
zone in early phases of the mixing (Fig. [3| panels b and c), 
it stays clearly separated later (Fig.[3j panel a). The blobs are 
always cooler and denser than the ambient matter by roughly 
0.4% and 0.05%, respectively (Fig.H]). Using the value of rel- 
ative density fluctuations Ap/p witnin the region where our 
mixing takes place we can derive an analytic estimate for the 
velocity of the dense blobs causing the finger-like structures, 
if we assume that the blobs arise from a dynamic instability 
( [Weiss et al.|2004) . With Ap/p ~ 5 x 10" 4 , inferred from our 



341 V 

simulation (Fig.|4]a), we find 



v ~ JgA— ~ 7 x 10 5 cm s" 



(1) 



where g ~ 10 8 cm s~ 2 is the gravitational acceleration at the 
base of the convection zone, and A ~ 10 7 cm the approxi- 
mate length of the fingers. The above estimate agrees well 
with the results of an analysis of our simulation, where the 
gas inside the fingers sinks at a bit smaller velocities rang- 
ing from ~ 1 x 10 5 cm s _1 up to 3.5 x 10 5 cm s" 1 . It confirms 
that the "buoyancy work" (~ v 2 ) is consistent with the accel- 
eration of gas seen in the simulation. While the initial back- 
ground state is dynamically stable (according to Ledoux and 
Schwarzschild), these flow properties suggests that the mix- 
ing takes place on a dyn amic timescale. This apparent contra- 
diction is addressed in jj |3.2| The development of the flow at 
t ~ 9 x 10 4 s is depicted in Figure^, d. 

Carbon Flash Convection. — A similar mixing process also 
appears in the carbon-flash model. A snapshot showing var- 
ious quantities in the lower half of the carbon burning shell 
convection zone and the underlying stable layer is presented 



Figs. p] and H which show the finger-like structures, the 
d in spherical coordinates (r and 6 in cm and degree) using 
a rectangular projection. This makes the finger-like structures appear longer 
than they really are. 



6 Note that in 
maps are displayed 
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Fig. 4. — Snapshots showing the base of the convection zone in the 2D model hefl.2d.3 at t = 63430 s. The panels (a,b,d) give the relative difference between 
the local and the horizontally averaged value at a given radius of the mean molecular weight //, the density p and the temperature T, respectively. The panels 
(c,e,f) display the deviation of angular velocity with respect to radius r dve/dr, radial velocity vr, and the Ledoux criterion (a positive value implies that the flow 
is Ledoux or dynamically unstable; see Eq.|2j, respectively. The vertical solid line marks the bottom of the convection zone (location of T max ), while the dashed 
line gives the location from where our mixing begins. Note that only a part of the computational domain is shown. 



in Figure Bb, e and can be seen to develop in a similar manner 
to the He flash model. 

In this model, however, the initial structure of the stable lay- 
ers is slightly different from that in the helium-flash model. 
When the model was mapped onto our hydrodynamic grid 
and allowed to adjust to hydrostatic equilibiurm, some narrow 
bands within the stable layers appeared which were unstable 
to convection. This was due merely to the differences between 
the stellar evolution code and hydrodynamics grid as well as 
the fact that the stable layer was only marginally stable. These 
narrow unstable regions quickly mixed, and stabilized, pro- 
ducing a slight stair stepping and shallow chemical disconti- 
nuities in the stable layer rather than a smooth profile as in the 
initial stellar model. Besides this initial transient, however, 
the evolution between this model and the helium flash model 
proceed in the same manner. 

Oxygen Shell Burning. — In contrast to the helium- and 
carbon-flash convection cases, the stable layer residing be- 
neath the oxygen shell burning zone does not undergo notice- 
able mixing over similar timescales. This is an interesting 
observation because the overall structure of the oxygen burn- 
ing shell is remarkably similar to the flash convection cases, 
except for the mean molecular weight gradient, which is pos- 
itive in the oxygen burning case. This fact strongly suggests 
that the molecular weight gradient is responsible, or at least, 
connected to the mixing seen in the shell convection mod- 
els. We also point out that the buoyancy frequency N 2 under 
convection shell is significantly higher than in the previous 
two cases indicating stronger dynamic stability (Figure [5]:, f), 
which could simply suppress the mixing. 



3.2. Linear Stability of the Mixing Layer 

As a starting point to analyze the mixing we consider the 
linear stability of the regions in question, the standard ap- 
proach to determine the mixing properties in ID stellar evo- 
lution codes. For adiabatic displacements, the condition for 
stability is often written in terms of the temperature and com- 
positions gradients, and reads 



AV < (ip/S)^ 



(2) 



where the stellar structure gradients are V = din T/d\n P, 
= dln/u/dlnP, the thermodynamic derivatives are 
V ad = (3 In 773 In JV, 6 = -{d mp/<9m T) P ^, <p = 
(d Inp/d In ia)p, s , and AV = V - V^. This stability condit ion is 
known as the Ledo ux criterion for convective stability (Kip- 
|penhahn & Weigert|1990| ) and it measures the degree to which 
the background is superadiabatic". In a homogeneous com- 
position medium the stability condition reads AV < 0, which 
is known as the Schwarzschild criterion and is equivalent to 
the statement that positive entropy gradients (ds/dr > 0) are 
stabilizing , negative entropy gradients are destabilizing, and 
zero entropy gradients result in neutral buoyancy upon dis- 
placement i n the linear regime when pressure equilibrium is 
maintained (Landau & Lifshitz 1959 ). The composition term 
in Eq. [2] accounts for the buoyancy due to mean molecular 
weight gradients in the medium. 

The stability condition (eq.(2)) can be identified with the 
region in the "nabla-plane" below the line described by AV = 
((f/S)V mu , as illustrated in Figure [6] Above this line the fluid 
is subject to dynamical instability (e.g., thermal convection, 
Rayleigh-Taylor instabilities), while below the line the fluid is 
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ditions (|Ulrich||1972[ |K ippenhahn et al. 1980; Grossman & 
Taam | 1 996[ |Charr7o nnel & Zahn 2007[ |Siess|2009[ |Stan cliffe 
et al.|2009D. 

Kippenhahn et al. (1980) first recognized that off-center 
helium and carbon ignition would indeed lead to thermoha- 
line mixing below the burning shells but the long (thermal) 
timescale estimated for the growth of this instability com- 
pared to the rapid evolution of the shell flash led to a dimin- 
ished interest in this context. In our shell flash simulations, 
the layers below the shell convection are indeed unstable ac- 
cording to the thermohaline condition (V < V a( i and < 0) 
but the mixing process operating in our simulations is in fact 
distinct from direct thermohaline mixing, since it operates on 
a significantly shorter timescale than the thermal timescale. 
Instead, we find that the regions which undergo mixing are 
forced by waves excited within the region by the adjacent tur- 
bulent convection zone. Therefore, what we have captured in 
our simulations is a wave induced, turbulent mixing process 
operating in a stably stratified layer. 



Fig. 6. — Temperature (AV = V - V a j) and composition gradient (0/5)V p 
plane. The region above the diagonal is dynamically unstable according to 
the Ledoux criterion (Eq.|2j. The lower right quadrant is dynamically sta- 
ble, while the "semiconvection" and "thermohaline" regions are dynamically 
stable, but unstable on thermal timescales. 

in a dynamically stable regime. But what about non-adiabatic 
effects? The regions in this plane labeled "semiconvection" 
and "thermohaline" are in fact unstable on thermal timescales 
due to secular instabilitie^] and their linear growth proper- 
ties have been analyzed for a sampling of astrophysical con- 

7 Secular instabilities are distinguished by growth timescales that are sig- 
nificantly longer than the dynamical timescale. 



3.3. Wave Driven Mixing 

As described in §3.1| and jj |3.2j the regions immediately be- 
low the connective shells are initially (linearly) dynamically 
stable. Furthermore, the absence of heat conduction in these 
simulations precludes these regions from undergoing ther- 
mally driven instabilities, such as semiconvection and ther- 
mohaline mixing. Therefore, the instabilities which develop 
in these regions are most naturally attributable to the jostling 
of the layers by turbulence in the adjacent convection zone. 
The nature of the jostling is predominantly in the form of in- 
ternal wave motions (i.e., gravity waves). 

As the stable layer is jostled by the neighboring convec- 
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tion, mixing is instigated which redistributes the entropy and 
the composition. This is illustrated in the time sequence pre- 
sented in Figure [7] for the helium shell flash model hefl.2d.3. 

There are several notable features in this figure. First, 
the convection zone remains separated from the underlying 
stable layer by a spike in buoyancy frequency N 2 (around 
r~ 4.7 x 10 8 cm in panels a and b) arising from the entropy 
jump between the convection zone and the underlying layer. 
This entropy jump is due to the energy deposition by the nu- 
clear burning in the shell. Note that the stabilizing layer is 
present despite the destabilizing \i -profile. In other words, the 
entropy profile is stable enough to compensate for the destabi- 
lizing yu-gradient effect. The spike in buoyancy also prevents 
overshoot mixing into the underlying stable layers. 

Another notable feature is the reduction in stability just be- 
low the buoyancy frequency spike which separates the region 
from the convection zone as time progresses. This reduction 
in stability (or decrease in N 2 ) is associated with a flattening 
of the entropy profile in the stable layer due to mixing in- 
duced in the region. The mixing that ensues is also associated 
with a modification of the composition profile, and by t ~ 10 4 
seconds a new step in the yu-profile has already formed near 
r ~ 4.5 x 10 8 cm. 

Finally, it is worth pointing out that in addition to the mix- 
ing inside the stable layer there is an overall deepening of the 
convection zone due to turbulent entrainment. By t ~ 9.5 x 10 4 
seconds the convective boundary has migrated inward a dis- 
tance of Ar ~ 2xl0 7 cm, fromr ~ 4.8xl0 8 cm tor ~ 4.6xl0 8 
cm. 

In the stable layers, the mixing is strongest in regions where 
the gradients are largest, as expected for a diffusive mixing 
process. Since compositional mixing at the molecular level is 



not included in these simulations, nor is heat conduction, this 
mixing results as a consequence of numerical diffusion. We 
stress here that the numerical diffusion can have a realistic 
counterpart in stars, but caution should be paid to the quanti- 
tative rate of mixing presented in this paper. We will discuss 
the connection between the simulated results (and numerical 
mixing) and the situation in real stars below in term s of wave 
enhanced diffusion as studied by, e.g., |Press] ( |1981| ). 



3.4. Layering and Interface Migration 

A striking phenomena seen in both the He- and C-flash con- 
vection models is the formation of homogeneous layers in the 
mixing region separated by steep gradient interfaces. These 
are apparent in both Figure [5] and Figure [7] Once these lay- 
ers have formed, a subsequent phase of evolution takes place 
which involves the migration of the interfaces separating these 
layers. This is illustrated in the series of snapshots compiled 
in Figure [8] which shows the time evolution of the buoyancy 
frequency N 2 in the mixing region for the carbon flash model. 
These N 2 profiles and their peaks are seen to evolve by a pro- 
cess of corresponding interface migration driven by turbulent 
entrainment. This entrainment is powered by the turbulent ki- 
netic energy deposited in the region by trapped gravity waves 
which are excited by the overlying convection zone, rather 
than thermally driven turbulence as in the case of a convec- 
tion zone. 

The formation of steps and interface migration is a well 
known phenomena in a variety of fluid dynamic systems, in- 
cluding those studied in geophysical systems, in laboratories, 
and in simulation. A very close analog to the type of inter- 
face migratio n that we see in our sim ulations is the situation 
modeled by Balmforth et al. ([1998, see their Figure 6 and 
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Fig. 8. — Profiles of the buoyancy frequency N 2 as a function of radial co- 
ordinate for different values of the time. The vertical axis is time and the 
profiles are shown spaced at time intervals of ~ 700 seconds. The inset figure 
at the top is showing the scale for N 2 corresponding to the last profile at time 
~ 53000 s. 

13, which show striking similarities to our Figure [8]). In their 
work, the evolution of buoyancy due to turbulence in a sta- 
bly stratified layer is modeled in an attempt to understand the 
origin of layer formation and the dynamics of interface mi- 
gration. This situation is a very close analog to the mixing 
process taking place in our simulations, except that the source 
of turbulence induced in the mixing region in our calculations 
is due to wave motions excited in the region. 

Although the quantitative rate of mixing in these regions 
remains to be identified through higher resolution calcula- 
tions and analysis, the formation of steps and their migration 
through the stable layer provide important clues about how 
one might treat this mixing in ID stellar evolut ion codes (e.g. 
IBalmforth et al.|1998l|ZTlitinkevich et al.|20Q7j ). 

4. DISCUSSION AND CONCLUSIONS 

We have identified a mixing process operating in our sim- 
ulations of off-center helium- and carbon-flash convection. 
The mixing results in the redistribution of composition and 
entropy within the stable layers that reside beneath the shell 
convection zones. Associated with the mixing process is the 
formation of "fingers" of high mean molecular weight ma- 
terial which traverse the mixing region. These "fingers" are 
essentially tracer of material trasport taking place across the 
stable layer (e.g., Figure [3}. 

We find that this mixing takes place against a background 



state which is initially dynamically stable to linear perturba- 
tions. And while the mixing region is indeed unstable to ther- 
mohaline mixing on a secular timescale, the mixing process 
that we observe operates on a dynamical timescale, thus ex- 
cluding thermohaline mixing as the root cause. This is re- 
inforced by the fact that heat conduction was not included 
in these calculations, primarily because its effects are too 
slow to operate over the time frames simulated. The thermal 
timescale for thermohaline mixing (Kippenhahn et al. 1980 ) 
estimated for our core helium flash model is of the order of 
10 5 years. It is neither expected that the timescale will drasti- 
cally change (decrease) with improved stellar models nor dur- 
ing th e flash as the thermal s tructures are in general very sim- 
ilar ( [Sweigart & Gross|1978l ). 

Our search for a counterpart to the mixing observed in 
our simulations among terrestrial experiments with qualita- 
tively similar conditions led us to phenomena like viscous or 
density fingering driven by a Rayleigh-Taylor instability ( |De| 
|Wit| |2004 , 2008) and buoyancy-driven instabilities resulting 
from differences in mass diffusion coefficients of chemically- 
reacting species sitting on top of each other ( Almarch a et al. 
2010). Properties of these phenomena do not fit the mixing in 
our simulations. 

As an alternative explanation we describe a scenario in 
which mixing is instigated by the finite amplitude jostling 
motions of the stable layer by the turbulence present in the 
adjacent convective shell. The perturbations induced in the 
stable layer by the adjacent convection are primarily in the 
form of waves (i.e., g-modes). The jostling of the stable layer 
results in a diffusive mixing process therein which leads to 
a redistribution of the entropy and composition on a dynam- 
ical timescales. This scenario is qualitativ ely s imilar to the 
enhanced mixing discussed by |Press|(|1981|) and|Press & Ry 



bicki ( 1981| ) who found that thermal and mass transport is not 
only a function of temperature and composition gradients, but 
also a function of the fluid shearing motions setup by the pres- 
ence of a passing wave field. 

Indeed, the layer from which the mixing begins coincides 
with layers of strongest shear (Fig .H] c). But the shear and 
the finite amplitude perturbations alone do not drive mixing 
and have to be preceded by the formation of sufficiently steep 
compositional gradients. In the case of our core helium and 
carbon flash models, such gradients are inherent in our initial 
stellar models at the bottom of convection shells due to off- 
center ignition. 

The biggest source of uncertainty inherent in this work is 
the quantitative rate of mixing. This uncertainty stems from 
the fact that the simulated mixing results from enhanced nu- 
merical diffusion, rather than from a resolved velocity field 
and molecular diffusion processes. Higher resolution simula- 
tions as well as analytic modeling should be pursued, perhaps 
in restricted models, in order to quantify the mixing rate due 
to this process. This is necessary for including it into a stellar 
evolution code. 

One of the pressing goals for stellar evolution is the con- 
struction of algorithms for ID codes that capture mixing pro- 
cesses such as that described in this paper, and which extend 
mixing due to hydrodynamic processes beyond the bound- 
aries of convective regions and into stably stratified layers. 
The need to treat turbulence and mixing in stable layers has 
already been recognized in the geophysic al community (e.g., 
|Fernando|1991[ |Zilitinkevich et al.|2007| ), and this work can 
serve to inform a related effort for stellar evolution. 
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